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Abstract 

In  this  paper,  a  fluid  model  of  the  incompressible  Stokes  equations  is  used  to  simulate 
the  motion  of  a  cell  boundary  separating  two  fluids  that  have  equal  or  different  viscosity 
and  non-linear  surface  tension.  Theoretically,  our  analysis  lead  to  decoupled  certain  jump 
conditions  that  are  used  in  constructing  the  numerical  algorithm.  Our  numerical  results 
agree  with  those  in  the  literature  and  include  some  new  results  that  may  can  be  used  to 
explain  the  cell  cleavage  process.  The  new  method  developed  here  preserves  the  area  very 
well  for  our  testing  problems.  Some  interesting  observations  are  obtained  with  different 
choices  of  the  parameters. 


1  Introduction 

In  the  literature,  the  Stokes  equations  are  widely  used  to  simulate  the  deformation  of  cells, 
see  [1,  2,  13]  and  many  others. 

In  reality,  cells  live  in  mediums.  Therefore  quite  often,  a  two  phase  model  is  reasonable. 
In  this  paper,  we  use  the  following  Stokes  equation  with  a  periodic  boundary  condition  to 
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study  the  motion  of  the  fluids  inside  and  outside  of  a  cell. 


Vp  =  pAu,  {x,y)  E  D  —  T, 

(1.1) 

V  •  u  =  0,  (x,y)  €  D, 

(1.2) 

Periodic  BC  for  u  and  p, 

(1.3) 

where  u  =  ( u ,  v)  is  the  velocity  vector,  p  is  the  pressure,  /i  is  the  viscosity,  which  we  assume 
to  have  constant  value  inside  and  outside  of  the  cell,  fl  is  a  bounded  domain,  T  is  the  cell 
boundary  which  is  closed  and  immersed  in  the  domain,  see  Fig.  1  for  an  illustration.  The 
periodic  boundary  condition  is  justified  since  the  environment  of  cells  does  not  vary  very 
much  and  widely  used  for  many  biology  problems,  see  [14,  15]  and  the  references  therein. 
Across  the  cell  boundary,  sometimes  is  called  the  interface,  there  are  the  following  jump 
conditions 


[p  —  2/r  (Vu  •  n,  Vv  ■  n)  •  n]  =  7 (s)k  (1.4) 

\H  (Vu  ■  n,  Vv  •  n)  ■  r]  =  l1-5) 

where  jump  [•]  is  defined  as  the  difference  of  the  limiting  value  from  outside  of  the  cell 
boundary  and  that  from  the  inside,  n  and  r  are  the  unit  normal  pointing  outwards,  and 
tangential  directions  of  the  cell  boundary  respectively,  7(5)  is  the  surface  tension  which  is 
the  concentration  of  tension  elements  and  is  modeled  as  a  separate  equation  in  [1,  2],  In  this 
paper,  it  does  not  have  to  be  a  constant  but  a  function  defined  on  the  surface  of  the  cell. 
The  surface  tension  is  the  main  driving  force  for  the  cell  deformation.  The  boundary  of  the 
cell  boundary  has  the  same  velocity  as  the  fluids  surrounding  it: 

dT 

dt 

Since  the  flow  is  viscous,  the  velocity  is  continuous  across  the  cell  boundary,  see  [4], 

Almost  all  the  numerical  simulations  for  cell  dynamics  with  a  free  boundary  using  the 
model  of  Stokes  equations  are  based  on  a  boundary  integral  method,  see  for  example,  [3, 
16,  17]  and  some  others.  Little  work  can  be  found  in  the  literature  using  the  finite  element 
method  or  the  finite  difference  method  to  solve  the  problems  because  of  the  almost  inhabited 
computational  cost  to  solve  the  problem  in  the  entire  domain  when  the  jump  in  the  viscosity 
is  large.  However,  when  the  viscosity  is  continuous,  the  problem  can  be  solved  easily  with 
the  finite  difference  method,  see  [8,  18]  and  others. 

In  this  paper,  we  use  a  finite  difference  method  based  on  a  fast  Poisson  solver  to  solve 
the  two-phase  moving  interface  problem.  One  of  advantages  of  this  approach  is  that  we  able 
to  obtain  the  velocity  and  pressure  in  the  entire  domain  without  sacrificing  the  speed  of  the 
simulation.  The  approach  described  here  is  not  tied  to  a  boundary  integral  equation  which 


(1.6) 
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Figure  1:  A  diagram  of  the  geometry  discussed  in  this  paper  and  it  is  used  for  the  derivation 
of  jump  relations. 

is  useful  to  construct  some  desired  force  along  the  cell  boundary  to  control  the  shape  of  the 
moving  interface  even  if  a  boundary  integral  equation  is  not  available  or  difficult  to  derive. 

Outside  and  inside  the  domain  excluding  the  boundary  of  the  cell,  the  equations  (1.1) 
can  be  reformulated  as  a  sequence  of  three  Poisson  problems,  one  for  each  variable  due  to 
the  incompressibility  condition.  For  piecewise  constant  viscosity,  applying  the  divergence 
operator  to  (1.1)  and  using  the  incompressibility  condition,  we  get 

Ap  —  0,  (x,y)efl  —  T,  periodic  BC  on  dfl  for  p,  (1.7) 

where  A  is  the  Laplacian  operator.  Since  the  right  hand  side  is  known,  this  is  a  Poisson 
problem  for  the  pressure.  Once  p  is  known,  the  equations  (1.1)  represents  two  independent 
Poisson  problems  for  u  and  v  respectively.  However  the  solutions  outside  and  inside  are  not 
independent.  A  numerical  algorithm  that  ignores  the  coupling  relations  usually  lower  the 
accuracy  especially  at  or  near  the  interface.  In  this  paper,  we  derive  the  jump  conditions  of 
the  solution  across  the  cell  boundary  and  use  them  to  construct  accurate  solution  method. 

The  rest  of  the  paper  is  organized  as  follows.  In  Section  2,  we  derive  the  jump  conditions 
for  (1 . 1 )- ( 1 .3) .  In  Section  3,  we  outline  our  numerical  algorithm  which  is  rather  straightfor¬ 
ward  if  we  know  the  jump  conditions.  Numerical  experiments  are  provided  and  analyzed  in 
Section  4.  Some  conclusions  will  be  given  in  Section  5. 

2  Decoupling  the  jump  conditions. 

Using  the  immersed  interface  method  (IIM),  [7,  9,  10,  11,  12],  we  can  solve  the  three  Poisson 
equation  efficiently  if  we  know  the  jump  conditions  in  the  solution  and  in  the  flux.  However, 
the  jump  conditions  given  in  (1.4)-(1.5)  are  coupled  and  we  can  not  use  them  directly.  The 
approach  used  in  [8]  for  Stokes  equation  can  not  be  applied  due  to  the  discontinuity  in  the 
viscosity.  Fortunately,  the  following  theorem  gives  decoupled  jump  conditions. 
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Theorem  1  Let  p,  u,  and,  v  be  the  solution  of  (1.1)-(1.3)  with  the  jump  conditions  (1-4)- 
( 1.5 ).  Then  the  following  equalities  hold 


w 

dl  (s) 
dr  T’ 

(2.8) 

b] 

=  7(s)  K, 

(2.9) 

bn] 

d27 

d^' 

(2.10) 

It  is  more  convenient  to  use  the  local  coordinates  to  prove  this  theorem.  Let  (X,Y)  be  a 
point  on  the  interface,  and  the  unit  normal  and  tangential  directions  are 

n  =  (cos  #,  sin#),  r  =  (—  sin#,  cos  9) 


respectively,  where  9  is  the  angle  between  the  outward  normal  direction  and  the  rr-axis,  see 
Fig.  2  for  an  illustration.  The  local  coordinates  at  (X,  Y)  then  is 


£  =  (x  —  X)  cos  9  +  [y  —  Y)  sin#, 

rj  =  —  (x  —  X)sm9+(y  —  Y)cos9. 

We  first  we  prove  the  following  lemma  which  is  needed  in  the  proof  of  Theorem  1. 


(2.11) 


V 


Figure  2:  Local  coordinates  and  geometry. 

Lemma  1  Let  u  be  the  solution  to  (l.l)-(l.S).  Then  the  following  equalities  are  true. 

[P  ( UXX  +  VXy)  ]  =  0,  [p  ( UXy  +  Vyy)  ]  =  0,  (2.12) 

[ pux  cos  9  +  vx  sin 9 ]  =  0,  [puy  cos  9  +  vy  sin#]  =  0.  (2.13) 

Proof:  In  the  inside  and  outside  of  the  domain  excluding  the  interface,  we  have 


ux  +  vy  =  0,  and  p  (ux  +  vy)  =  0, 

P  {'Uxx  T  ^xy)  =  0,  and  p  (uxy  +  v yy )  =  0. 
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Therefore  we  can  continuous  extend  the  definitions  above  to  everywhere  inside  the  domain 
and  we  have  the  equalities  in  (2.12)  right  way. 

To  prove  the  second  part,  we  enclose  the  interface  with  a  small  domain  fi(e),  see  Fig.  1  for 
an  illustration.  Let  <p(x,y)  E  C°°  be  a  two  dimensional  testing  function.  Since  y  ( uxx  +  vxy) 
and  ji  ( uxy  +  vyy)  are  continuous  inside  0(e)  after  the  extension,  we  have 


A4  ( uxx  +  vxy)  —  V  •  /r 


A4  (uXy  +  ^yy)  —  V  •  A4 


U1 


=  o, 


=  0. 


Here  we  have  used  the  fact  the  y  is  piecewise  constant  and 

,  \  d  (  du\  dfdv 

9\uxx  T  vxy)  —  )  +  7^3  (  A4-?:  I  —  0 


dx  \  dx )  dx  \  dyx 


(2.14) 

(2.15) 


has  continuous  extension  across  the  cell  boundary.  Multiply  (2.14)  by  any  testing  function 
4>(x,y)  E  C°°{ 0)  and  integrate  over  Of,  from  Green’s  theorem,  we  have 


lim 

e— >0 


(f)  yW  • 


u  x 
Vx 


dxdy 


/  cf)[y  (uxcos9  +  wxsin#)]  ds  +  lim  /  /  yVcf)' 

Jr  J  Jn £ 

J  cf)[y(ux  cos  9  +  vx  sin@)]ds, 


ux 

Vx 


dxdy 


where  we  have  used  the  fact  that  ux  and  vx  are  bounded.  Since  </;  is  arbitrary,  we  conclude 
that 


Similarly,  we  can  get 


[y  ( ux  cos  9  +  vx  sin  0)]  =  0. 


[y  ( uy  cos  9  +  vy  sin  0)]  =  0. 


(2.16) 


(2.17) 


That  completes  the  proof  of  the  lemma. 

Proof  of  Theorem  1.  Expressing  the  jump  condition  (1.5)  in  terms  of  the  local  coordi¬ 
nates,  we  get: 

(s') 

—  sintAI/ra^]  +  cos9[yv^\  +  cos9[yurj]  +  sin  9[yvv]  — - — — .  (2-18) 

Expressing  the  jump  condition  [y(ux  +  vy)\  —  0  in  terms  of  the  local  coordinates,  we  have: 

cos  9[yu^\  +  sin@[A«n^]  H —  sm9[yur]]  +  cos  9[yvv]  —  0  (2-19) 
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Under  the  local  coordinates,  the  two  equalities  in  (2.13)  are: 

cos2  6[/ju^\  +  sin  9  cos  0[/jv^]  —  s'm  6  cos  Ol/iUrj]  —  sin2  61/JVr)]  =  0,  (2.20) 

sin#cos 9\fJU(\  +  sin2  9\jj,v^\  +  cos2  Qlnu^]  +  sin#cos 9\jj,vv]  =  0.  (2-21) 

The  coefficient  matrix  of  (2.18)-(2.21)  is  non-singular  whose  inverse  evaluated  using  Maple 
is: 


(2.22) 


—  sin  9  0  2  cos2  9  —  1  2  sin  9  cos  9 

cos  9  0  2  sin  9  cos  9  —2  cos2  9+1 

0  —  sin  9  0  1 

0  cos#  —1  0 

Therefore  we  have  proved  that 

$7  chy 

[fm^\  =  sin#,  [fj,vz\  =  cos 9, 

which  is  (2.8)  in  Theorem  1.  The  proof  of  (2.9)  is  straightforward  if  we  substitute  the 
equalities  above  into  (1.4). 

From  the  system  of  equations  and  the  inverse  of  the  coefficient  matrix,  we  also  conclude 
that 


Imuv\  =  °>  =  0. 

Let  the  interface  T  be  £  =  x(rl)  in  the  neighborhood  of  (£,  rj)  =  (0, 0),  which  satisfies  x(0)  =  0, 
x'(0)  =  0.  Then  x,r(0)  —  K-  the  curvature  of  the  interface  at  (0,0).  Differentiating  [ur/]  =  0 
along  the  tangential  direction,  we  get 

[tiUrm]  =  [ 9vm ]  =  °-  (2-23) 

Differentiating  [//u^]  =  jyfr  along  the  tangential  direction,  we  get 

d2ry  d'y 

D*«cu]  =  q^2 T  +  «  frn-  (2-24) 

Expressing  [n(uxx  +  vxy)\  =  0  and  \/r(uxy  +  vyy)\  =  0  in  the  local  coordinates,  we  obtain 

[  9  {u&  cos2  $  —  2^^  cos  $  sin  $  +  ut]v  s™2  ^  cos  $  sin  $ 

+v^(cos2  9  —  sin2  9)  —  vm  cos  9  sin  9 )  ]  =0, 

[  n  (y,££  cos  9  sin  9  +  u^v  (cos2  9  —  sin2  9)  —  um  cos  9  sin  9 

+V££  sin2  9  +  2v^v  cos  9  sin  9  +  vvv  cos2  9 )  ]  =  0. 
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Multiplying  cos  9  to  the  first  equation  and  sin  9  to  the  second  equation  above,  after  some 
algebra,  we  get 


[  p  cos  9  +  V££  sin  9)  ]  =  [p  (u^v  sin  9  —  v^Tj  cos  9)  ] . 


(2.25) 


Now  we  are  ready  to  derive  the  jump  condition  for  ^  from  the  momentum  equation  in  the 
local  coordinates: 


[  Vp  •  n] 


[p  (u#  +  um)  ]  cos9  +  [p  (vtf  +  vvv)  ]  sin# 


[  p  (u£v  sin  9  —  V£V  cos  9)  ] 


d27 

dr'2 


sm 


<97 

~th 


527 

dr2 


2  9  —  k  ———  sin  9  cos  9  +  —4r  cos2  9  +  n  7-f  sin  9  cos  9 


<97 


d27 

(9  r2' 


This  completes  the  proof  of  Theorem  1. 


Remark  1  If  we  denote 


A  =  7,  h  =  ^  (2.26) 

as  the  normal  and  tangential  force  strength,  then  the  jum,p  conditions  proved  in  the  theorem 
are  the  exact  the  same  as  those  derived  in  [5,  8]  for  continuous  viscosity. 


3  Numerical  method 

For  simplicity,  we  assume  that  the  domain  is  a  rectangle  [a,  5]  x  [c,  d].  and  the  spatial 
spacing  is  h  =  (b  —  a)/M  =  {d  —  c)/N,  where  M  and  N  are  the  number  of  grid  points  in  the 
x  and  y  directions,  respectively.  A  standard  uniform  grid  is  used  but  the  scheme  discussed 
in  this  section  can  be  generated  to  other  grids  as  well  without  substantial  difficulty. 

We  use  a  marked  particle  approach  to  express  the  moving  interface  since  we  are  not 
particular  interested  in  topological  changes  at  this  stage.  The  advantages  using  the  particle 
approach  include:  (a)  taking  advantage  of  the  spline  interpolation  package  developed  in  [9] 
which  can  also  find  the  tangential  derivatives  easily  for  a  quantity  defined  on  the  interface, 
(b)  better  control  of  accuracy  and  area  preservation  since  we  have  control  of  numbers  of 
points  on  the  interface.  The  main  step  to  move  the  cell  boundary  from  time  level  tn  to  tn+1 
are  the  following: 

•  Evaluate  the  interface  information  such  as  the  tangential  and  normal  directions,  the 
jumps  and  their  derivatives  along  the  interface  using  the  period  spline  interpolation 
package  [9]. 
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•  Use  the  IIM  method  [7,  9,  10]  to  solve  the  Poisson  equation  for  the  pressure.  With 
the  IIM  method,  the  discrete  system  of  equations  are  standard  Laplacian  plus  some 
correction  terms  at  irregular  grid  points.  The  cost  is  just  a  little  more  than  one  call  to 
the  fast  Poisson  solver,  which  is  the  one  from  the  Fish-pack  in  our  method. 

•  Interpolate  the  pressure  to  get  px  and  py  including  those  grid  points  near  or  on  the 
interface.  The  standard  central  finite  difference  schemes  at  a  regular  grid  point  ( Xi,yj ) 
are 

„  _  Pi+lj  ~  Pi-lj  „  _  Pi, 3+ 1  ~  Pij~  1 

Px  ~  2  h  ’  py~  2  h 

At  an  irregular  grid  point,  we  use  the  one-sided  weighted  least  square  interpolation  [11] 
to  get  px  and  py.  For  example, 

Px{xi,yj)  =  ^PkiPki 
k,l 

where  the  summation  is  taken  for  those  points  (xk,yi)  that  are  in  the  same  side  of  the 
interface  as  (x,t,yj)  and 

{xk  -  Xi)2  +  (yi  -  y3  j2  <  5 2 

for  a  given  5  ~  Ch.  for  some  constant  C.  In  our  test,  C  is  usually  chosen  as  3.6/i.  Such 
interpolation  gives  smooth  error  distribution  and  better  accuracy.  Usually  the  system 
of  equations  for  the  interpolation  coefficients  are  under-determined  and  is  solved  using 
the  singular  value  decomposition  process.  Since  the  dimension  of  the  system  is  small, 
such  interpolation  does  not  add  much  extra  cost  to  the  entire  algorithm. 

•  Use  the  fast  IIM  method  [11]  to  solve  the  elliptic  equations  with  given  jump  conditions. 
Since  there  can  be  a  big  jump  in  the  viscosity,  such  a  solution  process  used  to  be  the 
most  expensive  part.  In  the  fast  IIM  method,  we  introduce  an  unknown  jump  \un] 
which  is  only  defined  along  those  marked  particle.  A  GMRES  iteration  is  used  to  solve 
the  unknown  \un] .  Each  iteration  of  GMRES  is  one  call  to  the  fast  Poisson  solver.  With 
the  fast  IIM  method,  the  number  of  iterations  is  independent  of  both  the  jump  in  the 
coefficients  and  the  mesh  size.  The  fast  IIM  method  takes  advantage  of  the  fast  Poisson 
solver  from  the  Fish-pack.  It  is  almost  impossible  to  simulate  the  problem  discussed 
here  with  large  jump  in  the  viscosity  without  the  fast  IIM  method.  The  fast  IIM  solver 
for  elliptic  interface  problems  with  piecewise  coefficients  is  available  through  the  public 
ftp  cite  at  North  Carolina  State  University1. 

•  Use  the  weighted  least  squares  interpolation  again  to  interpolate  the  velocity  field  at 
the  marked  points  to  get  the  velocity  U'Jj,  k  =  1,  2,  ■  ■  ■  n&. 


1ftp. ncsu.edu  under  the  directory  of  / pub / math / zhilin / Package. 
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Update  the  cell  boundary  using  the  computed  velocity 


X£+1  =  X£  +  ATU£,  k  =  1, 2,  ■  ■  ■  nb. 
The  time  step  for  the  explicit  method  is  chosen  as 

h2 


At  = 


max  |U^|  max{  /r+,  /r  } 


(3.27) 


(3.28) 


While  it  is  reasonable  to  use  an  explicit  method  for  two  dimensional  simulation,  an 
implicit  method  may  be  needed  for  three  dimensional  problems. 


4  Numerical  examples 

All  the  calculations  in  this  section  were  performed  at  North  Carolina  State  University  using 
Sun  workstations.  Most  simulations  are  done  within  hours  depending  on  the  number  of  grid 
points.  We  use  the  spline  interpolation  package  developed  in  [9]  to  express  the  interface. 

4.1  Example  1. 

We  first  test  the  classic  case  in  which  the  surface  tension  7  is  a  constant.  It  is  well  known 
that  the  cell  boundary  will  relax  to  a  circle,  its  equilibrium  state.  The  initial  cell  boundary 
is: 

r(0)  =  0.5  +  0.2  sin(50),  0  <  9  <  2tt,  (4.29) 

where  r  =  \J  x2  +  y'2.  The  computational  domain  is  [—1,  1]  x  [— 1,  1], 

We  test  our  code  for  different  viscosity  and  the  surface  tension.  The  numerical  experiments 
reveal  the  following  two  important  features.  In  Fig.  3,  we  plotted  the  evolution  of  the  cell 
boundary  with  a  80  by  80  grid  and  =  80  marked  points  on  the  cell  boundary. 

•  We  know  that  the  larger  the  surface  tension  is,  the  faster  the  boundary  restore  to  the 
equilibrium  state,  the  circle  with  the  same  area  as  the  original  cell.  In  Fig.  3,  we  plotted 
the  cell  boundary  at  different  time  with  different  surface  tension  and  the  ratio  /r_//r+. 
In  Fig.  3  (a)  and  (b),  the  surface  tension  is  7  =  0.01,  a  small  number,  and  it  takes  long 
time  for  the  cell  boundary  to  converge  to  or  close  to  its  equilibrium  state.  In  Fig.  3  (c) 
and  (d),  the  surface  tension  is  7  =  0.1,  we  see  the  cell  boundary  moves  much  faster 
than  the  first  case. 

•  Another  factor  that  affects  the  motion  of  the  cell  boundary  is  the  ratio  /r_//r+,  where 
/i+  and  ji~  are  the  viscosity  of  the  fluids  outside  and  inside  of  the  cell  boundary, 
respectively.  The  larger  this  ratio  is,  the  faster  of  the  motion.  In  Fig.  3  (a)  and  (c), 
the  ratio  is  one  and  we  see  slower  motion  compared  with  Fig.  3  (b)  and  (d)  even  if  the 
surface  tension  is  the  same. 
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The  algorithm  is  quite  accurate  in  the  sense  that  the  symmetry  and  the  area  are  both  pre¬ 
served  very  well.  In  all  cases,  the  relative  area  change  is  less  than  4%.  The  largest  area 
changes  are  from  first  a  few  steps  since  the  initial  cell  boundary  has  pretty  large  curvature. 
In  the  figure,  Earea  is  the  relative  error  in  the  area  at  that  particular  time. 

4.2  Example  2 

In  this  test,  we  try  to  use  a  non-linear  surface  tension  to  simulate  the  cell  cleavage  similar  to 
those  described  in  [1,  2,  3].  However,  we  are  dealing  with  two  phases  instead  of  one.  First 
we  test  the  case  that  is  similar  to  the  numerical  test  in  [3]  and  analyzed  in  [1 ,  2] .  The  initial 
cell  boundary  is  a  circle 

x2  +  y2  =  0.452. 

The  surface  tension  is  chosen  as 


r  =  0.01  +  e  5w(x(s)),  (4.30) 

where  x(s)  is  the  ^-coordinates  of  the  cell  boundary,  and  5w(x)  is  defined  as 

{- —  (1  +  cos  {-kx/2 w))  ,  if  \x\  <  2w, 

4wy  (4.31) 

0,  if  |a:|  >  2iu, 

the  discrete  delta  function  for  the  immersed  boundary  method  [14],  According  to  the  theory 
in  Greenspan  [1,  2]  and  verified  by  He  and  Dernbo  in  [3],  the  cell  will  bend  in  the  middle 
which  is  again  confirmed  by  our  numerical  example.  In  Fig.  4,  we  plotted  the  cell  boundary 
at  t  =  0  and  t  =  68.125  for  the  case  of  equal  viscosity,  and  t  =  18.2807  for  the  case  of  different 
viscosity  /r+  =  0.01  and  y~  =  1.  While  the  choice  of  the  non-linear  surface  tension  will  affect 
the  deformation  of  the  cell,  our  numerical  experiments  agree  with  the  conclusion  in  [3]  that 
the  cell  can  not  have  large  deformation  if  (4.30)  is  used.  In  fact,  for  two  dimensional  problem, 
the  boundary  will  stop  move  further  if  the  curvature  is  zero,  and  the  boundary  reaches  its 
equilibrium  state. 

In  order  to  have  a  cell  deform  further,  some  different  surface  tension  should  be  used  in 
the  model  (1.1)-(1.3).  We  use  the  following  surface  tension 

r  =  0.01  —  e  6w(x(s)).  (4.32) 

Note  that  this  is  similar  to  the  generalized  Gibbs-Thomson  condition,  see  [6]  and  the  refer¬ 
ences  therein.  The  initial  cell  boundary  is  the  circle  r  =  0.45  with  a  perturbation: 

r{9)  =  0.45  +  0.15  cos  20,  0  <  6  <  2tt.  (4.33) 

This  can  be  interpreted  as  the  follows:  In  some  small  part  of  the  cell  boundary,  there  is 
sudden  change  in  the  concentration  of  the  material  due  to,  for  example,  chemical  reaction. 
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Figure  3:  The  cell  boundary  at  different  time  computed  using  a  80  by  80  grid  and  a  constant 
surface  tension  7  =  0.01  in  (a)  and  (b),  and  7  =  0.1  in  (c)  and  (d).  While  the  boundary 
in  all  cases  approaches  to  the  equilibrium  state,  a  circle,  the  time  scale  is  very  different. 
The  boundary  moves  faster  as  7  or  the  ration  /r_//r+  gets  larger,  (a)  The  result  with  equal 
viscosity  [i~  =  /r+  =  1,  and  7  =  0.01  at  time  t  =  0,  t  —  39.4531,  t  =  109.7656,  and 
t  =  163.28.  The  shape  is  far  away  from  the  equilibrium  state  even  with  t  =  163.  (b)  The 
result  with  different  viscosity  /r“  =  1,  n+  =  0.01,  and  7  =  0.01  at  time  t  —  0,  t  =  14.5128, 
t  =  37.9503,  and  t  =  85.1378.  The  shape  is  very  close  to  the  equilibrium  state  at  t  =  85. 
So  the  motion  is  much  faster  than  the  case  with  equal  viscosity,  (c)  The  result  with  equal 
viscosity  n~  =  /r+  =  1,  and  7  =  0.1  at  time  t  —  0,  t  —  10.6074,  t  =  30.1222,  and  t  =  77.0128. 
The  motion  is  faster  than  the  case  with  small  surface  tension  but  slower  than  the  case  with 
larger  ratio  of  n~  / /r+ .  (d)  The  result  with  different  viscosity  n~  =  1,  /r+  =  0.01,  and  7  =  0.1 
at  t  =  0,  t  =  1.0761,  t  =  2.9288,  and  t  =  14.6779.  The  motion  is  faster  than  any  other  three 
cases  in  this  figure. 
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(a) 


(b) 

M=160,  Nb=160,  H=(1,1)  M=160,  Nb=160,  p+=0.01,  g  =1 


Figure  4:  The  cell  boundary  at  different  time.  The  dotted  plot  is  the  initial  cell  boundary 
which  is  the  circle  r  =  0.45.  The  surface  tension  is  given  in  (4.30).  (a)  The  solid  line  is  the 
cell  boundary  at  t  =  68.125  in  the  case  of  equal  viscosity  /r+  =  n~  =  1.  The  relative  area 
change  at  t  =  68.125  is  0.64%.  (b)  The  solid  line  is  the  cell  boundary  at  t  =  18.2807  is  the 
case  of  different  viscosity  ji+  =  0.01  and  ji~  =  1.  The  relative  area  change  at  t  =  18.2807  is 
0.59%.  Again  the  motion  is  faster  when  the  ratio  /r_//r+  is  larger. 


(a) 


Evolution  of  the  front,  M=160,  Nb=160,  g=(1,1) 


(b) 

Blow-up  of  tips,  M=160,  Nb=160,  p=(1 ,1) 


Figure  5:  The  evolution  of  the  cell  boundary  for  for  Example  3  with  equal  viscosity  ji+  = 
=  1  at  different  time,  (b)  Blow-up  of  the  upper  tip.  The  relative  area  changes  at  different 
time  are  also  listed  in  the  figure. 
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As  the  result,  the  cell  will  gradually  deform  as  observed  in  the  experiments.  In  Fig.  5, 
we  plotted  the  evolution  process  of  the  cell  boundary  with  equal  viscosity  /i+  —  /i~  —  1, 
w  —  5/z,  and  e  =  h,  where  h  —  (b  —  a)  /M  —  (d  —  c) /N.  is  the  spatial  step  size.  We  see  that 
cell  is  dividing  and  eventually  will  break  up.  The  simulation  is  done  with  a  160  by  160  grid 
and  160  marked  points  on  the  cell  boundary.  Fig.  5  (b)  is  the  blow-up  plot  of  the  upper  tip. 
The  area  changes  which  is  also  marked  on  the  figure,  are  relatively  small,  less  than  3.2%,  for 
such  a  problem  with  the  tips  deepening  so  much  in  both  directions.  In  Fig.  6,  we  plotted  the 
simulation  with  different  viscosity  //+  =  0.01  and  //,  =  1,  we  see  similar  behavior  as  the  case 

of  the  equal  viscosity  but  faster  motion  as  we  discussed  for  Example  1. 


(a) 


Evolution  of  the  front,  M=160,  Nb=160,  g=(0.01,1) 


(b) 

Blow-up  of  tips,  M=160,  Nb=160,  g=(0.01 ,1) 


Figure  6:  The  evolution  of  the  cell  boundary  for  for  Example  2  with  different  viscosity 
//+  =  0.01  and  //,  =1  at  different  time,  (b)  Blow-up  of  the  upper  tip.  The  relative  area 

changes  at  different  time  are  also  listed  in  the  figure.  The  motion  is  faster  than  the  case  of 
the  equal  viscosity. 


5  Conclusion. 

We  have  developed  a  finite  difference  method  for  Stokes  equations  for  incompressible  two- 
phase  flow  that  models  cell  deformation  with  discontinuous  viscosity  and  non-linear  surface 
tension.  The  decoupled  jump  conditions  for  the  pressure  and  the  velocity  across  the  cell 
boundary  are  crucial  for  the  success  of  the  new-developed  algorithm.  With  similar  choice  of 
the  surface  tension,  our  numerical  results  agree  with  the  results  and  analysis  in  the  literature. 
Our  numerical  results  also  showed  that  with  suitable  choice  of  surface  tension,  a  cell  can 
deform  and  break-up.  While  it  is  known  that  motion  of  a  cell  boundary  depends  on  the 
surface  tension,  our  numerical  experiments  also  showed  that  the  motion  depends  on  the 
difference  of  the  viscosity  as  well. 
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Future  work  we  intend  to  study  including:  (a)  design  implicit  scheme  so  that  we  can  use 
larger  time  steps,  (b)  Three  dimensional  and  axisymmetric  simulations,  (c)  Optimal  choice 
of  the  surface  tension  to  match  certain  patterns  of  cell  deformations. 
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